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O ■ ABSTRACT 

w 

^ ' A 2.5-dimensional particle-in-cell code is used to investigate the propaga- 

tion of a large-amplitude, superluminal, nearly transverse electromagnetic (TEM) 
^ I wave in a relativistically streaming electron-positron plasma with and without 

I a shock. In the freestreaming, unshocked case, the analytic TEM dispersion 

relation is verified, and the streaming is shown to stabilize the wave against 
parametric instabilities. In the confined, shocked case, the wave induces strong, 
coherent particle oscillations, heats the plasma, and modifies the shock density 
profile via ponderomotive effects. The wave decays over > 10^ skin depths; the 
decay length scale depends primarily on the ratio between the wave frequency 
and the effective plasma frequency, and on the wave amplitude. The results 
are applied to the termination shock of the Crab pulsar wind, where the de- 
cay length-scale (> 0.05") might be comparable to the thickness of filamentary, 
variable substructure observed in the optical and X-ray wisps and knots. 



Subject headings: pulsars: general — pulsars: individual (Crab Pulsar) — stars: 
neutron — shock waves — plasmas — outfiows 
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1. INTRODUCTION 

A rotation-powered pulsar emits most of its spin-down luminosity as a wind of rel- 
ativistic particles, mainly electrons (e~) and positrons (e"'') produced in pair cascades in 
charge-starved regions of the magnetosphere (Arons 2002). Observations of pulsar winds 
confined by supernova remnants (Hester et al. 2002) or the interstellar medium (Chatterjee 
& Cordes 2004; Gaensler et al. 2004) reveal that the wind terminates in a synchrotron- 
emitting shock with a cylindrically symmetric ('crossbow') morphology, implying that the 
outflow is coUimated geometrically along the rotation axis, while the energy flux is concen- 
trated in the equatorial plane (Komissarov & Lyubarsky 2004). The wind is launched from 
the light cylinder with most of its luminosity transported as Poynting flux (Coroniti 1990), 
and its enthalpy flux is negligible beyond the light cylinder tl due to adiabatic cooling. On 
the other hand, models of the shock kinematics, the postshock flow, the synchrotron spec- 
trum (Kennel & Coroniti 1984b), and the variability of shock substructure (e.g. wisps and 
knots in the Crab) (Spitkovsky & Arons 2004), strongly suggest that the ratio a of Poynting 
flux to kinetic-energy flux is small at the shock, with a ~ 10~^ in several objects (Kennel & 
Coroniti 1984a; Gaensler et al. 2002). 

In recent years, theoretical attention has focused on the radial structure of the wind, in 
order to explain how a decreases from the magnetosphere {a » 1) to the shock (cr <^ 1): 
the (7-paradox. The dissipation of Poynting flux seems to be related to the decay of the 
large-amplitude wave launched into the wind at the pulsar spin frequency. One possible 
cause is that ideal magnetohydro dynamics (MHD) breaks down beyond a critical radius rc, 
in the sense that there are insufficient charge carriers to screen out the rest-frame electric 
field, and the electron inertia term in the relativistic Ohm's law becomes dominant (Melatos 
& Melrose 1996; Melatos 1998; Gedalin et al. 2001; Kuijpers 2001; Melatos 2002). In this 
respect, particular attention has been paid to the entropy wave carrying the alternating 
('striped') magnetic field in the wind (Bogovalov 1999; Coroniti 1990) and its dissipation 
by magnetic reconnection at the corrugated current sheet separating the magnetic stripes 
(Lyubarsky & Kirk 2001), possibly emitting observable high-energy pulses (Arons 1979; Kirk 
et al. 2002; Kirk 2004; Skjajraasen 2004; Petri & Kirk 2005). The dissipation occurs upstream 
from the termination shock, if the reconnection rate and particle fiux are large enough (Kirk 
& Skjaeraasen 2003), or else in the shock itself (Lyubarsky 2003). 

In this Letter, we assess an alternative scenario in which, beyond Tc, the plasma up- 
stream from the termination shock carries a superluminal transverse electromagnetic (TEM) 
wave, which modulates the wind and the shock (Melatos 2002). The wave is taken to be 
monochromatic, with a frequency given by the pulsar spin rate; an assumption generally 
used also in MHD wind models. Coherent particle acceleration in the TEM wave typically 
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leads to (7 <S 1; a self-consistent, WKB wave model predicts a ~ 10~^ at the termination 

shock given the measured particle and energy fluxes of the Crab pulsar (Melatos & Melrose 
1996; Melatos 1998). We use the 2.5-dimensional, relativistic, electromagnetic particle-in-cell 
(PIC) code XOOPIC (Verboncoeur et al. 1995) to investigate the free propagation of the 
wave upstream from the shock (§2), plasma heating in the shock (§3), and the fate of the 
wave inside and downstream from the shock (§3). We apply our results to the Crab pulsar 
wind in §4. 



2. FREELY PROPAGATING TEM WAVES 

A nonlinear TEM wave can propagate even if its frequency uo is less than the effective 
plasma frequency ^J j e^vri'^^ ^ where is the number density in the laboratory 
frame, m is the electron mass, and is the Lorentz factor, because the wave forces charges 
to oscillate relativistically, increasing their effective mass. In a cold plasma, the plane wave 
dispersion relation is (Akhiezer & Polovin 1956; Kaw & Dawson 1970; Kennel & Pellat 1976; 
Melatos & Melrose 1996) 

= 1 (1) 



where r] — ck/u \s the index of refraction, ujp — -^/rv^e^ / e^m'^d, la is the bulk Lorentz factor, 
and is a dimensionless amplitude defined by = eEo/mccu, where Eq is the physical 
electric field amplitude. For E > 1, the charges are accelerated from nonrelativistic to 
relativistic speeds within a single wave period. Note that a nonlinear electromagnetic wave 
is strictly transverse only if it is circularly polarised. In this case, the streaming speed cf3d 
is a constant of the motion and can be chosen arbitrarily. As a result, a oc E'^'-jcT^ varies 
independently with E and 7^^, whereas in a linearly polarised nonlinear wave, 7^; (and hence 
(t) is determined uniquely by E. In the small-amplitude limit, n^^'^d and E decouple. 

In order to verify (1) with XOOPIC, we launch a cold beam from the left-hand 
edge {x = 0) of the simulation box. Two orthogonal, phased dipole antennas are placed at 
X — 0, continuously emitting a circularly polarised TEM wave with frequency uj — ujq which 
ramps up to a constant amplitude E over a rise time tr < 0.2tc, where tc is the time for 
the beam to cross the box. The right-hand edge of the box is transparent both to particles 
and waves (thanks to a wave-absorbing algorithm). Once initial transients disappear, we 
Fourier transform the transverse electric field component Ez{x,t) and measure 7] for the 
Fourier component of largest amplitude. For linear waves with 10"^ < E < 10^^, we verify 
the linearized version of (1) to approximately 0.5 per cent in vacuo, 1 per cent in a uniform. 
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cold stationary plasma, and 2 per cent for a relativistically drifting, uniform, cold plasma. 
In the evanescent regime u; < ujpy/2, the wave penetrates a few plasma skin depths c/up. 

When a TEM wave with E ^ 1 is launched into a nonstreaming plasma, its pondero- 
motive force excavates a density cavity around the antenna. If the energy density of the 
wave field, eoEg, exceeds n^mc^, the cavity expands at relativistic speeds. To verify (1) for 
a nonstreaming plasma, one would need to set up the wave field and particle distribution 
self-consistently at all x before starting the simulation (which we defer to a future paper), 
and even then, one would expect the wave to be disrupted by parametric instabilities (Max 
& Perkins 1972; Sweeney & Stewart 1975; Ashour-Abdalla et al. 1981; Leboeuf et al. 1982). 
However, if a nonlinear TEM wave is launched into a streaming plasma (7^ ^ 1), it prop- 
agates freely — without excavating the plasma and without being disrupted by parametric 
instabilites — in the regime u > ijjp\/2/(l + E'^y/^. In the parameter ranges 1 < E < 10^ 
and 1 < 7d < 3.9 x 10^, the value of rj we compute agrees with (1) to 5 per cent. We stress, 
however, that more work is needed to rigorously verify that the simulated fields and particle 
momentum distributions are compatible with the assumption of an infinite, plane TEM wave 
given by (1). 



3. SHOCKED TEM WAVE 

In certain respects, the above PIC experiment with freely propagating TEM waves over- 
laps with previous works (Ashour-Abdalla et al. 1981; Leboeuf et al. 1982). The main focus 
of the current work is a different scenario: A nonlinear TEM wave which encounters an ultra- 
relativistic shock formed by a confining medium, as in the case of a pulsar wind termination 
shock. In order to simulate such a shock, we repeat the experiment described in §2 but add 
a particle-reflecting, wave-transmitting magnetic wall in the region 95 < 5; = xutp-^/c < 100, 
where Up^ = ujp{x = 0). We trace 0(10^ — 10^) particles for 0(10^ — 10^) timesteps on a grid 
where each cell has a size Ax — Ay < O.lcTr/cUp-f, and the timestep is given by At — OAAx/c. 
The grid has 128 cells and periodic open boundaries in the y-direction, and from 256 to 1600 
cells in the x-direction. Below, the antenna (TEM wave) frequency is denoted by coq. We 
observe that the after reflected off the magnetic wall, the counterstreaming particles trigger 
a Weibel instability, which causes the beam to fllament in the y direction. The initial stage 
of the instability is marked by exponentially growing oscillations (along the y axis) of the 
magnetic fleld component Bz, with characteristic wavelength c/up. A shock forms after a 
time t given by u!p-^{t — t^.^^) ~ 40 — 70, where t^.^^ is the instant when the front of the 
beam reflects off the magnetic wall. As there is no dc magnetic fleld in our simulations, the 
cyclotron instability seen by others (Gallant et al. 1992; Hoshino et al. 1992; Spitkovsky & 
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Arons 2004) is absent. 

In the discussion below, we average all quantities over y, smoothing out the filaments. 
We fix the bulk Lorentz factor at x = to be 70 = 3870. The value of a at the antenna, a{x — 
0), is denoted by Eq, and we adopt parameters such that E(x = 0) = 5.4 x 10^(ujp^/uJo)y/Eo. 
Downstream from the antenna, the field couples to the particles and quickly reduces E until, 
at i; ~ 25, a 'fully developed' TEM wave forms, in which E(x) ~ u±{x)/c, where u±{x) is 
the local ensemble average of (uy"^ + u.^'^Y^'^ (uy and are components of the 4-vclocity). 
We refer to the regions 10 < £ < 50, 50 < £ < 80, and a; > 80 as the shock precursor, the 
shock interior and the downstream medium, respectively. 

3.1. Plasma Heating 

In the shock precursor, there is a strong, coherent, position-dependent interaction be- 
tween the wave and the particles. The lower panels of Fig. 1 show the e+ number density 
in (£, Uz)-spa.ce at two different positions. At low £, the particle motion is almost phase- 
coherent, but electrostatic and electromagnetic fluctuations gradually increase the thermal 
spread with x until it exceeds the wave-induced quiver motion. The conversion of Poynting 
fiux and kinetic-energy fiux into thermal motions is evident from the upper panels of Fig. 
1, which show the frequency spectrum of the electric field component Ey. As x increases, 
the thermal background (broad sidebands) grows at the expense of the antenna-driven TEM 
mode (narrow lines at a; = ±<^o)- 

The positron energy distribution 7(7) is shown at two positions in Fig. 2, one in the 
shock precursor and the other in the shock interior. Interestingly, for a wave with Eq > 1, we 
find (7) (X 70 (Eo -I- 1) at both positions (indeed, throughout the shock), although the form 
of 7(7) changes. In the precursor (8 < x < 16), the distribution is narrow and drops quickly 
at high energies; fitting the tail with a power law, 7(7) ^ 7^", wc obtain 15 < a < 35 
for the range 2 < Sq < 24. This does not necessarily imply nonthermal acceleration, but 
rather the heating of an initially cold distribution. In the shock interior (50 ^ x ^ 58), /(7) 
turns into an anisotropic Maxwellian (see §3.3). The stronger the TEM wave, the hotter the 
downstream plasma, with (7) determined by the total upstream energy fiux. 

The anisotropy of the particle distribution in the shock interior depends on the po- 
larization of the wave; for example, a linearly polarized wave with E^ ^ 1 and Ey = 
predominantly energizes Uz- A detailed study of this issue is outside the scope of this paper, 
which focuses on circularly polarized TEM waves. 
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3.2. Fate of the Wave 

As the wave propagates away from the antenna and accelerates the particles into oscil- 
lations, a quickly drops below Sq. For x > 25, E changes little until it encounters the shock, 
beyond which the wave continues to decay. The decay rate scales with the skin depth, c/up^. 
Figure 3a shows the TEM wave amphtude as a function of position and Eq. For Eq = 0.4, 
we find E{x = 2S)/E{x = 0) = 0.05 and E{x = 4:0)/E{x = 0) < 0.01. For Eq = 3.6 (14.0), 
the corresponding ratios are 0.20(0.25) and 0.07(0.14), respectively. In the case Eq = 14.0, 
E remains roughly flat for 40 < i; < 70, before dropping away for x > 70. The drop is 
accompanied by a temperature decline, possibly because the plasma near the magnetic wall 
is not yet thermalized. For Eq — 6, a similar drops occurs at x > 57. 

Taken together with similar simulations spanning the parameter range 1 < ouo/cup^ 10, 
Fig. 3a leads to an important result: the length scale Ld over which E drops by two orders 
of magnitude, scales proportionally to the skin depth as 

Ld ^ Kc/iVp^, (2) 

where 10^ ^ K < lO'^ depends on uiq/uJpY, i-e., near (far from) from the cutoff given by (1) the 
decay is faster (slower). This is consistent with previous results for an unshocked, relativistic 
EM wave (Leboeuf et al. 1982). The value of K also increases slowly with the amplitude 
(for 0.2 < So < 20 we obtain 50 < < 150 with cuo/u^pf — 2.0). For the parameter regime 
<-^o/<-^pt which we cannot achieve in our simulations (since it would require a very large 

number of timesteps to both allow the shock to form and resolve the wave oscillations) , the 
value of K might exceed 10^. Over distances much greater than L^, the wave may enter the 
linear regime, but we are unable to distinguish the wave from the noise at this level in our 
simulations. (Note that if Ld is redefined as the length scale over which E drops by a factor 
of e, K is reduced by a factor of 4.6). 

As seen from Fig. 3b, as the wave decays the tranverse momentum density e±^{x) — 
n{x)u±{x)/c, where n{x) — n'^{x) +n~{x), increases with x. The largest-amplitude case 
(So = 14.0) is an exception, with de±/dx < near x = and the magnetic wall; this may 
be due to ponderomotive forces (§3.3), imperfect absorption of the wave by the simulation 
boundaries, or other boundary effects. 

The wave decay might be caused by a parametric instability which sets in when the 
bulk speed of the flow is reduced in the shock interior. Further simulations are needed to 
obtain enough statistics to rigorously settle the microphysical aspects of the decay; this is 
deferred to a forthcoming work. 
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3.3. Density Profile 



Our code tracks all three vector components of the electromagnetic fields, but only 
supports electrostatic (Langmuir) fiuctuations in two directions. Particles therefore dif- 
fuse slower in than in and Uy. In this respect, our simulations are essentially two- 
dimensional, with an effective adiabatic index F ~ 3/2 (not 4/3, as in isotropic 3D). This 
explains the form of the density profile n{x) in Fig. 4, plotted t — t^Q^ + SbcOp^'^ for three 
different wave amplitudes. In the near-hydrodynamic case Eq — (not shown), the density 
jump is n2/ni ~ 2.8. where rii {n2) denotes the upstream (downstream) density, consistent 
with the Rankine-Hugoniot prediction n2/ni = F/(r — 1) ~ 3 for F = 1.5. Conversely, 
given the observed density jump for Eq — 0, we infer F = 1.56, implying some slow diffusion 
in Uz- Figure 4 also suggests that ^2/^1 increases with Sq. This trend can be ascribed to 
the ponderomotive force executed by the decaying TEM wave, and/or by the longitudinal 
electric field (x E induced by parametric decay to the Langmuir mode. As the TEM wave 
decays, we have a ponderomotive force proportional to dE'^/dx; hence ^2/^1 increases with 
Eq. Interestingly, we see little or no evidence that the shock speed changes with Eq. 



CRAB PULSAR WIND 



For a steady-state pulsar wind, conservation of mass and energy give 
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In (3) and (4), ujq is given by the pulsar spin rate, the total particle injection rate is = 
10^^s"^7V38, the solid angle filled by the wind is Q-^, and we define /j, — L/Nmc^, where L 
is the pulsar spin-down luminosity. For canonical Crab parameters (5 x 10^ A* ^ 5 x 10^; 
77 ~ 1, rtyj/An ~ 0.1 - 1, a{rs) ~ 3 x 10"^ and r^/r^ ~ 2.7 x lO'^) (Triimper & Becker 
1998; Weisskopf et al. 2000; Kirk & Skjaeraasen 2003; Spitkovsky & Arons 2004), we take 
oJpt = ^p{rs) and find 1 < E < 10, 9 x 10^ < cuo/up^ < 3 x 10^ and 9 x 10^ < L^/r^ < 3 x 10« 
at the termination shock (using 10^ % K < 10^). 

In a self-consistent model of the TEM wave, a can be expressed in terms of the other 
quantities in (1), (3), and (4) according to Eq. (97) of Melatos & Melrose (1996), viz. 

1 cu^il- rf) 



,7/(1-77') 47,2a;4 



(5) 
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This allows us to solve self-consistently ior rj,a,E, and (Vp-^ given just N and 7^. Adopting the 
current best estimates of these quantities, N^s = 10°'^ and 7^ = 10^, as inferred from an ion- 
cyclotron model of the variability of the Crab wisps (Spitkovsky & Arons 2004), together with 
~ 47r, we solve (1) and (3)-(5) to obtain cr ~ 9.0 x 10~^, E ~ 2.5, L^/r^ ~ W -W, and 
1 — 7/ ~ 4.8 X 10^^^ (the value for t] reflects the fact that the displacement current dominates 
the conduction current for r > Tc). 

Evidently, our simulations explore artificially small values of uJo/oOp'^ (for best-guess 
values, Wo/^^pf ~ 9 x 10^ at the Crab shock), and artificially large values of a and E (e.g., 
in Fig. 4, Sq = 3.6 gives a ~ 0.01 and ~ 4 x 10^ at x ~ 50) due to numerical limitations. 
Nevertheless, we consistently find ^ lO^c/wp-j- for 0.1 Sq < 100, i.e., over three decades 
in Sq. The fate of the TEM wave thus seems clear: It decays beyond r^, just hke the entropy 
wave (Lyubarsky 2003), but the dissipation mechanism is not reconnection. The decay scale 
exceeds 10^ r^. 

Observationally, might be resolved. It subtends 0.05" — 0.5" if the Crab is at distance 
of 2 kpc, using 10^ < K < 10^ and the best estimates for the wind parameters. This range 
includes the angular thickness of filamentary, variable substructure in the equatorial wisps 
of the Crab pulsar wind, observed by the Hubble Space Telescope and the Chandra X-ray 
Observatory (Hester et al. 2002), but at present the connection is speculative. Since 
increases with uo/up^ in the parameter regime near cutoff, and the Crab shock value of uj/up^ 
cannot be achieved in our simulations, we can only conclusively give a lower limit (lOOc/wpi) 
for L(i; i.e., the possibility that > lO^c/Up-^ (> 0.5") cannot be excluded. A theoretical 
upper limit for can be estimated by noting that, ultimately, adiabatic cooling of the 
diverging, subrelativistic flow beyond the Crab termination shock enforces a cutoff, and 
parametric instabilities are bound to set in, yielding < 0{rs)- 

The TEM wave generates a latitude-dependent, nonthermal particle distribution inside 
the shock which is anisotropic near the equatorial plane, where the wave is (approximately) 
linearly polarized, and isotropic at high latitudes, where the polarization is (approximately) 
circular. In principle, future X-ray and/or optical polarization measurements can test this 
prediction. 

The simulations reported here do not include a dc magnetic field, as exists for an oblique 
rotator. Such a field affects the propagation and stabifity of the TEM wave (Asseo et al. 
1980). For example, a subluminal, circularly polarized TEM wave with a dc magnetic field 
in the x-direction (radially) is stable in the inner part of the wind but unstable to three- wave 
parametric decay in the outer wind (Melatos 1998, 2002). These important issues will be 
addressed in a forthcoming paper. 
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Fig. 1. — The shock precursor. Top row: Frequency spectrum |£'y(ci;)| at (a) x = 14 and 
(6) X = 29, obtained during the time interval 56 < (^pf(t — ^^-ef) ^ 80- The hues at u = 
±lvq are the nonhnear TEM wave. As x increases, the hues decay and broaden while the 
thermal background spectrum grows. The lines at ^ iO.Gtuo are due to electrostatic 
plasma oscillations. Bottom row: Positron density in {x, Uz)-spa.ce for (c) 10 < x < 20 and 
(d) 20 < x < 30. In (c), the particles are strongly phase-coherent with the wave (which 
has a wavelength of fa 3.3x), although thermal broadening can be seen at x > 11. In 
(d), the thermal width is larger than the quiver amplitude, but a periodic modulation is 
still evident (in the absence of the wave, one gets a straight beam centered on = 0, 
gradually broadening with x). The parameters used are Eq — 2.0, 70 — 3870, a;o/<^pt = 1-97, 
t = + and E{x = 0) = 2.0 x IQ^VEq. 
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Fig. 2. — The energy distribution 7(7) of positrons for 8 < x < 16 (a) and 50 < x < 58 (b), 
normalized such that J d^f{^) — n'^. The labels denote Eq; the other parameters are as in 
Fig. 1. 




Fig. 3. — (a) Fourier amplitude of the TEM wave at a; = a;o as a function of x, normalized 
to its value at x = 0. The solid line represents free propagation (no shock). The other 
curves arc for shocked TEM waves. Arrows indicate upper limits, and labels denote Eq- 
(b) Transverse momentum density e±{x) = n{x)u±{x)/c, in arbitrary units. As the wave 
decays, e± increases due to plasma heating. Each curve represents a shocked TEM wave. 
Both panels are snapshots at i ~ t^^^ + lOOcUp^^^. Other parameters are as in Fig. 1. 
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Fig. 4. — The aggregate density profile n{x) — + for a TEM wave shock. The 

regions 10 < x < 50, 50 < x < 80, and x > 80 are the shock precursor, the shock interior and 

the downstream medium, respectively (see §2). The drop at a; > 95 is due to the magnetic 
wall. The solid, dashed-dotted, and dashed lines are for Eg = 24, Eq = 2.0, and Eq = 0.4, 
respectively. Other parameters are as in Fig. 1. 



